---
title: "Mt. Vesuvius Lives"
subtitle: "INFO 526 - Summer 2025 - Final Project"
author:
- name: "Seismic Minds - Weston Scott"
affiliations:
- name: "School of Information, University of Arizona"
description: "Exploratory Analysis of Seismic Activity: Mt. Vesuvius"
format:
html:
code-tools: true
code-overflow: wrap
code-line-numbers: true
embed-resources: true
editor: visual
code-annotations: hover
execute:
warning: false
echo: false
render-on-save: true
---
```{r library_setup, include=FALSE}
if (!require("pacman"))
install.packages("pacman")
# Use pacman::p_load to install and load CRAN packages
pacman::p_load(
countdown,
# datasauRus,
dplyr,
elevatr,
gganimate,
ggplot2,
ggthemes,
ggspatial,
gifski,
glue,
grid,
gridExtra,
here,
hms,
knitr,
lubridate,
metR,
openintro,
readr,
rnaturalearth,
rnaturalearthdata,
rnaturalearthhires,
terra,
scales,
sf,
tidyverse,
viridis
)
# Handle GitHub package separately
if (!require("dsbox")) {
# Install devtools if not present
if (!require("devtools"))
install.packages("devtools")
devtools::install_github("tidyverse/dsbox")
library(dsbox)
}
ggplot2::theme_set(ggplot2::theme_light(base_size = 12))
options(width = 65)
knitr::opts_chunk$set(
fig.width = 7, # 7" width
fig.asp = 0.618, # the golden ratio
fig.retina = 3, # dpi multiplier for displaying HTML output on retina
fig.align = "center", # center align figures
dpi = 300 # higher dpi, sharper image
)
save_figs <- TRUE
```
## Data Wrangling
```{r}
vesuvius_data <- read.csv(here("data","vesuvius.csv"))
vesuvius_data_filtered <- vesuvius_data |>
drop_na(
latitude,
longitude,
depth_km,
duration_magnitude_md
) |>
mutate(
depth_abs = abs(depth_km),
depth_boundary_1km = case_when(
depth_km < 1 ~ "< 1 km",
depth_km < 2 ~ "< 2 km",
depth_km < 3 ~ "< 3 km",
depth_km < 4 ~ "< 4 km",
depth_km < 5 ~ "< 5 km",
TRUE ~ "> 5 km"),
depth_boundary_1km = factor(
depth_boundary_1km,
levels = c("< 1 km",
"< 2 km",
"< 3 km",
"< 4 km",
"< 5 km",
"> 5 km")),
depth_boundary_0.5km = case_when(
depth_km < 0.5 ~ "< 0.5 km",
depth_km < 1 ~ "0.5 ~ 1 km",
depth_km < 1.5 ~ "1 ~ 1.5 km",
depth_km < 2 ~ "1.5 ~ 2 km",
depth_km < 2.5 ~ "2 ~ 2.5 km",
depth_km < 3 ~ "2.5 ~ 3 km",
depth_km < 3.5 ~ "3 ~ 3.5 km",
depth_km < 4 ~ "3.5 ~ 4 km",
depth_km < 4.5 ~ "4 ~ 4.5 km",
depth_km < 5 ~ "4.5 ~ 5 km",
TRUE ~ "> 5 km"),
depth_boundary_0.5km = factor(
depth_boundary_0.5km,
levels = rev(c("< 0.5 km",
"0.5 ~ 1 km",
"1 ~ 1.5 km",
"1.5 ~ 2 km",
"2 ~ 2.5 km",
"2.5 ~ 3 km",
"3 ~ 3.5 km",
"3.5 ~ 4 km",
"4 ~ 4.5 km",
"4.5 ~ 5 km",
"> 5 km")))
) |>
filter(
!is.na(depth_boundary_0.5km)
) |>
mutate(
depth_rescaled = rescale(depth_km, to = c(1, 11)),
mag_rescaled = rescale(duration_magnitude_md, to = c(0.4, 0.9))
) |>
arrange(
depth_boundary_0.5km,
desc(depth_rescaled)
) |> glimpse()
```
## Display Basic Stats
```{r}
max_depth <- max(vesuvius_data$depth_km, na.rm = TRUE)
min_depth <- min(vesuvius_data$depth_km, na.rm = TRUE)
print(paste("Maximum depth (km):", max_depth))
print(paste("Minimum depth (km):", min_depth))
max_mag <- max(vesuvius_data$duration_magnitude_md, na.rm = TRUE)
min_mag <- min(vesuvius_data$duration_magnitude_md, na.rm = TRUE)
print(paste("Maximum Magnitude:", max_mag))
print(paste("Minimum Magnitude:", min_mag))
```
## Make Initial Scatter Plot
```{r}
original <-
ggplot(
vesuvius_data_filtered,
aes(x = longitude,
y = latitude)
) +
geom_point(
aes(color = duration_magnitude_md,
alpha = depth_abs),
size = 3
) +
scale_color_viridis_c(
option = "magma",
name = "Magnitude"
) +
scale_alpha_continuous(
name = "Depth (km)",
range = c(0.75, 0.25)
) +
scale_x_continuous(
limits = c(14.34, 14.48),
expand = c(0, 0)
) +
scale_y_continuous(
limits = c(40.78, 40.88),
expand = c(0, 0)
) +
labs(
title = "Seismic Events at Mount Vesuvius",
subtitle = "Point Color by Magnitude\nPoint Transparency by Depth",
x = "Longitude",
y = "Latitude",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot",
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.direction = "horizontal",
legend.box = "vertical", # stack legends on top of each other
legend.box.just = "left",
legend.spacing.x = unit(0.3, "cm"), # tighter spacing
legend.spacing.y = unit(0.3, "cm"), # tighter spacing
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.background = element_rect(fill = "transparent", colour = NA),
legend.box.background = element_rect(fill = "transparent", colour = NA)
) +
guides(
color = guide_colorbar( # Magnitude
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 1),
alpha = guide_legend(
direction = "horizontal",
title.position = "top",
label.position = "bottom",
override.aes = list(alpha = c(0.75, 0.5, 0.25)),
keywidth = unit(0.4, "cm"),
keyheight = unit(0.4, "cm"),
order = 2)
)
original
if (save_figs) {
ggsave("images/original_vesuvius_plot.png", original)
}
```
## Generate Histogram for Depth Analysis
```{r}
histogram <-
ggplot(
vesuvius_data_filtered,
aes(x = duration_magnitude_md)
) +
geom_histogram(
binwidth = 0.2,
fill = 'red',
color = "white"
) +
facet_wrap(
~ depth_boundary_0.5km,
scales = "free_y"
) +
labs(
title = "Histogram of Duration Magnitude",
subtitle = "Distribution Facets by Depth",
x = "Magnitude (md)",
y = "Count",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot"
)
histogram
if (save_figs) {
ggsave("images/histogram_vesuvius_plot.png", histogram)
}
```
## Capture Elevation Data for Contour Mapping
```{r}
# Define Vesuvius bounding box (in WGS84 lon/lat)
vesuvius_bbox <- matrix(c(14.33, 40.78, 14.48, 40.89), ncol = 2, byrow = TRUE)
# Convert bbox to sf polygon
vesuvius_sf <- st_as_sf(
st_sfc(st_polygon(list(rbind(
c(vesuvius_bbox[1,1], vesuvius_bbox[1,2]),
c(vesuvius_bbox[2,1], vesuvius_bbox[1,2]),
c(vesuvius_bbox[2,1], vesuvius_bbox[2,2]),
c(vesuvius_bbox[1,1], vesuvius_bbox[2,2]),
c(vesuvius_bbox[1,1], vesuvius_bbox[1,2])
))), crs = 4326))
# Get elevation data from AWS Terrain Tiles (zoom 10)
elev_raster <- get_elev_raster(
locations = vesuvius_sf,
z = 10,
clip = "locations")
# Convert to dataframe for ggplot
elev_df <- as.data.frame(
elev_raster,
xy = TRUE)
colnames(elev_df) <- c("x", "y", "elevation")
glimpse(elev_df)
```
## Create Static Display of Depth Bins
```{r}
static_plot <- ggplot() +
metR::geom_contour_fill(
data = elev_df,
aes(x = x,
y = y,
z = elevation/1000,
group = NULL),
bins = 10
) +
scale_fill_gradientn(
colors = gray.colors(10, start = 1.0, end = 0.3),
name = "Elevation (km)",
limits = c(-0.111, 1.246), # Match elevation range from elev_df (km)
breaks = seq(0, 1.2, by = 0.3)
) +
geom_point(
data = vesuvius_data_filtered,
aes(x = longitude,
y = latitude,
color = duration_magnitude_md,
size = depth_rescaled),
alpha = 0.75
) +
scale_size_continuous(
name = "Depth (km)",
breaks = c(2.5, 5, 7.5),
range = c(1, 11) # adjust to control point size
) +
scale_color_viridis_c(
option = "magma",
name = "Magnitude (md)"
) +
scale_x_continuous(
limits = c(14.34, 14.48),
expand = c(0, 0)
) +
scale_y_continuous(
limits = c(40.78, 40.88),
expand = c(0, 0)
) +
labs(
title = "Seismic Events at Mount Vesuvius",
subtitle = "Point Color by Magnitude\nPoint Size by Depth",
x = "Longitude",
y = "Latitude",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot",
axis.title = element_text(size = 14),
axis.text = element_text(size = 12),
legend.direction = "horizontal",
legend.box = "vertical", # stack legends on top of each other
legend.box.just = "left",
legend.spacing.x = unit(0.3, "cm"), # tighter spacing
legend.spacing.y = unit(0.3, "cm"), # tighter spacing
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.background = element_rect(fill = "transparent", colour = NA),
legend.box.background = element_rect(fill = "transparent", colour = NA)
) +
guides(
fill = guide_colorbar( # Elevation
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 1),
color = guide_colorbar( # Magnitude
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 2),
size = guide_legend( # Depth
direction = "horizontal",
title.position = "top",
label.position = "bottom",
override.aes = list(alpha = 0.7),
keywidth = unit(0.15, "cm"),
keyheight = unit(0.4, "cm"),
order = 3),
alpha = "none"
)
static_plot
if (save_figs) {
ggsave("images/static_vesuvius_plot.png", static_plot)
}
```
```{r include=FALSE}
ggplot2::theme_set(ggplot2::theme_light(base_size = 14))
options(width = 65)
knitr::opts_chunk$set(
fig.width = 7, # 7" width
fig.asp = 0.618, # the golden ratio
fig.retina = 2, # dpi multiplier for displaying HTML output on retina
fig.align = "center", # center align figures
dpi = 150 # higher dpi, sharper image
)
```
## Create Animated Display of Depth Bins
```{r}
vesuvius_anim_plot <- ggplot() +
geom_point(
data = vesuvius_data_filtered,
aes(x = longitude,
y = latitude,
color = duration_magnitude_md,
size = depth_rescaled),
alpha = 0.75
) +
scale_x_continuous(
limits = c(14.34, 14.48),
expand = c(0, 0)
) +
scale_y_continuous(
limits = c(40.78, 40.88),
expand = c(0, 0)
) +
scale_size_continuous(
name = "Depth (km)",
breaks = c(2.5, 5, 7.5),
range = c(1, 11)
) +
scale_color_viridis_c(
option = "magma",
name = "Magnitude (md)"
) +
labs(
title = "Seismic Events at Mount Vesuvius",
subtitle = "Point Color by Magnitude\nPoint Size by Depth: {closest_state}",
x = "Longitude",
y = "Latitude",
caption = "Source: Seismic data for Mount Vesuvius\nTidyTuesday 2025-05-13\nGraphic: Weston Scott"
) +
theme(
plot.title.position = "plot",
plot.title = element_text(margin = margin(b = 23),
size = 16),
axis.text = element_text(size = 12),
legend.direction = "horizontal",
legend.box = "vertical", # stack legends on top of each other
legend.box.just = "left",
legend.spacing.x = unit(0.3, "cm"), # tighter spacing
legend.spacing.y = unit(0.3, "cm"), # tighter spacing
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
plot.caption = element_text(margin = margin(t = 30), size = 10, hjust = 1),
legend.background = element_rect(fill = "transparent", colour = NA),
legend.box.background = element_rect(fill = "transparent", colour = NA)
) +
guides(
color = guide_colorbar(
direction = "horizontal",
title.position = "top",
barwidth = unit(2.8, "cm"),
barheight = unit(0.4, "cm"),
order = 1),
size = guide_legend(
direction = "horizontal",
title.position = "top",
label.position = "bottom",
override.aes = list(alpha = 0.7),
keywidth = unit(0.15, "cm"),
keyheight = unit(0.4, "cm"),
order = 2),
alpha = "none"
) +
transition_states(
depth_boundary_0.5km,
transition_length = 3,
state_length = 2
) +
shadow_mark(
past = TRUE,
future = FALSE,
alpha = 0.2
) +
enter_fade() +
exit_fade()
gganimate::animate(vesuvius_anim_plot)
if (save_figs) {
gganimate::anim_save("images/vesuvius_animation.gif", vesuvius_anim_plot)
}
```
